// -----------------------------------------------------------------------------
// M A R I T I M E   R E S E A R C H   I N S T I T U T E   N E T H E R L A N D S
// -----------------------------------------------------------------------------
// Copyright (C) December 15, 2009 MARIN - All rights reserved - http://www.marin.nl
// -----------------------------------------------------------------------------
// Program  :
// Module   :
// File     : PulseManager.cpp
// -----------------------------------------------------------------------------

#include <new>

#include "defines.h"
#ifndef _WINDOWS
#include "mermaid.h"
#endif
#include "PulseManager.h"
#include "PulseFactory.h"
#include "Pulse.h"
#include "Geometry.h"
#include "assertplus.h"

#define COLUMNS  (unsigned)(Nx1 + 3)
#define ROWS     (unsigned)(Nx2 + 3)

const float oneSixteenth = 1.0 / 16;
const float oneEighth    = 1.0 / 8;
const float oneFourth    = 1.0 / 4;

// -----------------------------------------------------------------------------
// function PulseManager::PulseManager
// -----------------------------------------------------------------------------
/// constructor
/** detailed description */
// -----------------------------------------------------------------------------
PulseManager::PulseManager(int Nx1,
                           int Nx2,
                           REAL dx1,
                           REAL dx2,
                           REAL pstart,
                           bool smoothing)
: m_Nx1(Nx1)
, m_Nx2(Nx2)
, m_dx1(dx1)
, m_dx2(dx2)
, m_start(pstart)
, m_smoothing(smoothing)
, m_pulses()
, m_pulseField(COLUMNS, ROWS)
, m_pulseFieldOld(COLUMNS, ROWS)
{
    m_pulseField = 0;
    m_pulseFieldOld = 0;
}


// -----------------------------------------------------------------------------
// function PulseManager::~PulseManager()
// -----------------------------------------------------------------------------
/// destructor
/** detailed description */
// -----------------------------------------------------------------------------
PulseManager::~PulseManager()
{
    unsigned int n = (unsigned int)m_pulses.size();
    for (unsigned int i = 0; i < n; ++i)
    {
        delete m_pulses[i];
    }
}


void PulseManager::addPulse(Geometry* a_pGeometry)
{
    Pulse *pPulse = PulseFactory::create(a_pGeometry, m_dx1, m_dx2);
    m_pulses.push_back(pPulse);
}


unsigned int PulseManager::getNumberOfPulses() const
{
    return (unsigned int)m_pulses.size();
}

Pulse* PulseManager::getPulse(int a_nr) const
{
    ASSERT(a_nr < (int)getNumberOfPulses(), "no such pulse");
    return m_pulses[a_nr];
}

void PulseManager::setSmoothing(bool on)
{
    m_smoothing = on;
}

void PulseManager::lower(int a_nr, REAL a_currenttime)
{
    getPulse(a_nr)->setStartTime(a_currenttime);
}

void PulseManager::computePulsePattern(REAL a_currenttime, REAL dt,
                                       const Array2D<REAL>* pNetHeight,
                                       Array2D<REAL>** ppPressurePattern)
{
    bool timesmoothing = dt > 0.000001;
    if (timesmoothing)
    {
        m_pulseFieldOld = m_pulseField;
    }
    REAL g = getGravity();
    m_pulseField = 0;
    int n1 = m_pulseField.dim1();
    int n2 = m_pulseField.dim2();
    unsigned int n = (unsigned int)m_pulses.size();
    for (unsigned int ii = 0; ii < n; ++ii)
    {
        Array2D<REAL> pulse;
        int x;
        int y;
        m_pulses[ii]->compute(&x, &y, &pulse);
        // pulse contains depth data of ship in current position and attitude
        int ni = (int)(pulse.dim1());
        int nj = (int)(pulse.dim2());
        for (int i = 0; i < ni; ++i)
        {
            for (int j = 0; j < nj; ++j)
            {
                if (pulse[i][j] < -0.001) // there is a certain depth
                {
                    int xx = x + i + 1; // Klopman arrays have ...
                    int yy = y + j + 1; // ... offset 1 !
                    if (xx >= 0 && xx < n1 && yy >= 0 && yy < n2)
                    {
                        REAL wave = pNetHeight ? ((*pNetHeight)[xx][yy]) : 0;
                        REAL netPulse = wave - pulse[i][j];
                        REAL pslow = 0.0;
                        REAL starttime = m_pulses[ii]->getStartTime();
                        if (starttime > 0.0)
                        {
                            pslow = pSlow(a_currenttime - starttime, m_start);
                        }
                        REAL pressurePulse = pslow * g * netPulse;
                        // pressure pulse data must be positive!
                        if (pressurePulse > 0)
                        {
                            m_pulseField[xx][yy] += pressurePulse;
                        }
                    }
                }
            }
        }
    }
    if (timesmoothing)
    {
        REAL F = 2 * PI * dt / sqrt(2 * PI * std::min(m_dx1, m_dx2) / g);
	    for (int i = 0; i < n1; ++i)
	    {
		    for (int j = 0; j < n2; ++j)
		    {
                m_pulseField[i][j] = F * m_pulseField[i][j] + (1 - F) * m_pulseFieldOld[i][j];
            }
        }
    }
    if (m_smoothing) // spatial smoothing
    {
        /* smoothing pattern:
         * 1/16   1/8   1/16
         * 1/8    1/4   1/8
         * /16    1/8   1/16 */
        
        // copy pulsefield to temporary store
        Array2D<REAL> tmpPulseField = m_pulseField;
        m_pulseField = 0;
        
        // fill pulsefield again from temporary store while smoothing
        for (int i = 1; i < n1 - 1; ++i)
        {
            for (int j = 1; j < n2 - 1; ++j)
            {
                m_pulseField[i][j] = oneSixteenth *(tmpPulseField[i - 1][j - 1] + 
                                                    tmpPulseField[i - 1][j + 1] +
                                                    tmpPulseField[i + 1][j - 1] +
                                                    tmpPulseField[i + 1][j + 1]) +
                                     oneEighth    *(tmpPulseField[i - 1][j] +
                                                    tmpPulseField[i + 1][j] +
                                                    tmpPulseField[i][j - 1] +
                                                    tmpPulseField[i][j + 1]) +
                                     oneFourth    * tmpPulseField[i][j];
            }
        }
    }
    *ppPressurePattern = &m_pulseField;
}

void PulseManager::getGridProperties(int* pNx1, int* pNx2, REAL* pdx1, REAL* pdx2)
{
    ASSERT(pNx1, "pNx1=0");
    ASSERT(pNx2, "pNx2=0");
    ASSERT(pdx1, "pdx1=0");
    ASSERT(pdx2, "pdx2=0");
    *pNx1 = m_Nx1;
    *pNx2 = m_Nx2;
    *pdx1 = m_dx1;
    *pdx2 = m_dx2;
}
